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The Apollonian packings (APs) are fractals that result from a space-filling procedure with 
spheres. We discuss the finite size effects for finite intervals s £ [smin, Smax] between the largest 
and the smallest sizes of the filling spheres. We derive a simple analytical generalization of the 
scale-free laws, which allows a quantitative study of such physical fractals. To test our result, a new 
efficient space-filling algorithm has been developed which generates random APs of spheres with a 
finite range of diameters: the correct asymptotic limit Smin/smax — >■ and the known APs' fractal 
dimensions are recovered and an excellent agreement with the generalized analytic laws is proved 
within the overall ranges of sizes. 
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1. Introduction 

The problem of finding the circle inscribed into the in- 
terstices between mutually tangential circles and tangent 
to them (a so-called osculatory packing), historically at- 
tributed to Apollonius of Perga, was solved by Descartes 
and independently rediscovered various times 1 . Leib- 
niz pointed out the possibility of obtaining a peculiar 
kind of packing by iterating the procedure of inserting 
such kissing circles, whose size decreases as the insert- 
ing procedure goes on [2|- By starting from an initial 
configuration and by recursively filling the space with 
the osculatory packing down to arbitrarily small diam- 
eters, the Apollonian packing (AP) is constructed. In 
Fig. [T] two examples of AP are presented for Euclidean 
dimension d — 2. The structure of an AP is self-similar 
because it is repeated on different scales of observation. 
In general, the self-similarity can be exact or statistical 
and leads to a fractal. The main quantity characteriz- 
ing a fractal structure is the fractal dimension, df, which 
is a (Lipschitz) invariant of the set descending from the 
Hausdorff-Besicovitch (HB) measure definition and pos- 
sibly differs from the Euclidean dimension d O H] • 
The osculatory packing construction has been extended 
to d > 2 and the study of systems of polydisperse spheres 
is made in the attempt of understanding how their frac- 
tality affects some macroscopic observables and, by re- 
verse, how their formation mechanism influences their 
fractality. In particular, in d = 3, the AP has been pro- 
posed as a model for dense granular systems [5]. It has 
been also used in describing stress yielding properties 
of materials, for example in concretes as well as in 
the study of seismic gaps or geological faults [7^. More- 
over the scale-free properties of the AP are of particular 
interest in the context of complex networks [H |^ . Re- 
cently AP of non-spherical objects have also been studied 

While the application and the characterization of the AP 
has been widely studied, there is no exact theoretical pre- 
diction for the value of d/ in 2 and 3 dimensions and 



various techniques have been developed to build AP and 
numerically evaluate their df. For example, it has been 
calculated for the plane tiled with circles obtained by the 
circular inversion method [12|. In c? = 3, a generalized 
inversion algorithm has been adopted J^j. 
So far we have discussed the deterministic AP but it has 
been proved that its fractal nature emerges also when a 
random sequence of space-filling insertions is pursued. In 
this generalized model [T^ [T3] , called random Apollonian 
packing (RAP), the circles are inserted one at time with 
center positions randomly chosen; after the insertion, the 
diameter is simply inflated until it touches a previously 
inserted circle. Extended models make use of simultane- 
ously inflating circles. All these models, deterministic or 
not, are shown to have universal features belonging to a 
broader class of models called "packing-limited growth" 
|15j . The RAP relies on the fact that all the osculatory 
packings in a certain dimensionality must have the same 
fractal dimension [16]. Several routes to RAP has been 
devised; we suggest Refs. [TTJ [TSl [T7] as a short review. 




FIG. 1: Apollonian packings for d — 2: a) classical example 
as in Ref. [18]; b) an example of random packing from the 
algorithm (smax/i ~ 1) as explained later in the text. 

Both deterministic and random procedures, when stud- 
ied numerically, are affected by finite-size effects. In the 
context of AP (as in its random counterpart), starting 
from a given configuration of equal spheres of diameter 
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Smax, the impossibility of having packings with arbitrar- 
ily small sizes means that the asymptotic limit 



^ 



(1) 



(where df is defined) will never be reached. 
It would be interesting to systematically characterize 
the deviation from such a limit. In this paper, we 
propose a simple solution to this problem and we 
find that, taking finite size effects into account, the 
fractal dimension remains well defined also for finite 
ranges s S [smin, Smax]- The paper is organized as 
follows. In the next two paragraphs, we introduce the 
basic ideas and observables and we provide an analytic 
solution for the problem of evaluating df for a finite 
size range. In the following sections, we will introduce 
a new algorithm to construct osculatory random AP in 
arbitrary dimension with the possibility of tuning the 
size ratio defined in eq. ([I]). After having verified its 
correct asymptotic behaviour, we shall finally use this 
novel algorithm to test the proposed finite size correction. 



2. Apollonian packings and physical fractals 

In the rest of the paper, we refer to the hard spheres 
as the set of non- overlapping geometrical objects (also 
for dimensions different from d = 3) and the size s then 
refers to the diameter. Given a packing of equal hard 
spheres which occupy a volume fraction (/>, the porosity 
£. = 1 — (p can be decreased by filling the free interstices 
with smaller spheres; by iterating the insertion procedure 
of smaller spheres, the final packing results to be a frac- 
tal as Sniin gocs to 0, which corresponds to the limiting 
value 0=1. At equal number N of total inserted objects, 
the Apollonian packing (AP) is known to be the densest 
of these packings [6]. Despite the universality of d/ for 
all osculatory packings [TB], i.e. for all configurations 
in which any interstice is filled with the largest possi- 
ble hard sphere, there is not yet an analytical expression 
of the AP's fractal dimension. The fractal dimensions 
df — 1.305684 and df = 2.473946 have been numerically 
calculated for the AP in d = 2 and d = 3 respectively 

[HIIS]. 

In general the relation d— l<df<dis valid, as con- 
firmed by numerical simulations 19 and the AP scale- 
free nature is expressed by the size distribution 



n(s) oc s-^'^^+i) , 



(2) 



a power-law defined for all the positive diameters s [H]. 
Fractals are implicitly understood to be scale-free struc- 
tures and eq. ([2| is a consequence of that. However, when 
a finite range of sizes exists, the fractal is called a phys- 
ical fractal [20j|2Tj. In this case, the distribution of s is 
limited to an interval [smim •Smax]7 where Sniin and s^ax 
are the smallest and the largest diameters of the pack- 
ing, respectively. Various practical methods have been 



introduced to calculate the fractal dimension: the box- 
counting, for example, is applied in the case of physical 
fractals [22l [23] and to direct measurements of physical 
systems |24) . 

We use two typical observables, the inverse cumulative 
distribution (number of spheres with size greater or equal 
to Smin) and the porosity, which under the asymptotic 
condition ([T]) are respectively given [12] by 



A^(Smin) = ^ 1 OC sj^ , 
£(Smin) = 1 - 



d—df 



(3a) 
(3b) 



These measures can be used to estimate the df for 
packings of polydisperse hard spheres: for different oc- 
currences of Smin, the fractal dimension is computable as 
smaller Smin are considered in measuring the values ( 3a 



and ( 3b I [6 . All the different methods for evaluating the 



d^ of packings rely on the fact that the fractal dimension 
definition works only approaching the asymptotic condi- 
tion 

Q- This means that, for finite ranges of diameters, 
relevant deviations exist starting from Smin/smax 1/5, 
as explicitly highlighted in Ref. [5]. 



2.1 Finite size deviation 

The main idea is that the geometrical building rule itself, 
with its iterativity, defines the fractal-like behavior, while 
the finite interval of sizes influences only the quantity of 
objects used for the building. The self-similarity of a 
set then implies the existence of a set of similarities in 
the generation of the fractal, which is beyond the mere 
agreement of the value of the fractal dimension [25] . 
This idea can be quantitatively rendered. For a random 
AP where the osculatory packing constrain is respected 
at each insertion (building iterative rule) we expect that 
the value for df will remain the same. Due to the fractal 
nature of AP and RAP the power-law ([2| should hold 
and we make the ansatz that the finite size correction is 
completely accounted for the proportionality constant of 
the distribution. In the finite case, we rewrite n(s) as 
nf{s): 



-idf + l) 



(4) 



The recursive fractal construction is accounted by the 
power-law and we now calculate the corrections to 
eqs. (3a) and (3b I through the use of /(smin, Smax)- 
The fraction of space occupied by Ng^^^^ spheres of max- 
imum size is by definition 



V 



(5) 
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where V is the total volume and the curvature factor is 

tor even a 
for odd d 



(J _ '^^'^''^ _ I 2^ (d/2)\ 
" ~ T{d/2)d 1 1 



2('i-i)/2 rf!! 

The total number of spheres N and volume fraction (p 
can be expressed using the 0— th and d— th moments of 
the distribution Q as: 

N = Ns„^^^^ + nf{s')ds' , (6a) 

= + (Cd/V) ^ '^^^\f{s')sMs' . (6b) 

A natural condition arises assuming that all the space 
shall be occupied as the filling procedure continues; this 
means that, for any positive value of Smaxj it is 



lim ^(j) — 1 



(7) 



Using eq. Q into eq. (6b I with the limit ([T]), we obtain 
the distribution's proportionality constant as 



— {d- d/)Sma: 



(8) 



It is important to notice that it does not depend on Sn 



Inserting this result into eqs. (6a) and (6b) finally gives: 



N 



id - df) 



d—df 



, (9a) 



(9b) 



Note that eqs. (3a) and (3b) are the particular asymp- 



totic cases of their more general expressions (9a I and 
(9b), as expected. These results remain valid in the 



limit V oo. Despite the simple hypothesis made, now 
we have an explicit expression for the observables which 
allows to evaluate the deviation from the ideal case ([!]). 

2.2 Filling algorithm 

In order to test the previous results, AP have been gener- 
ated with the help of a new numerical algorithm (which 
works in any Euclidean dimension d > 1). The developed 
algorithm has the same basic behavior as the random 
Apollonian packing (RAP) mechanism, where the filling 
process starts with an initial population of hard-spheres 
of a specified diameter (smax) and proceeds with new 
spheres added one at a time into the unoccupied space; 
randomly fixing the center of any new sphere, its size is 
determined by extending the diameter until it touches 
its closest sphere. But the RAP is not expected to build 
osculatory packings at each insertion, as we demand in 
order to check the theoretical previsions. In the algo- 



rithm we propose, random movements are additionally 
performed by the inserted spheres, in order to enhance 
the filled space accepting only the displacements that al- 
low its diameter to grow (possibly up to s,„ax)- 
No distribution of sizes for the spheres to insert is a pri- 
ori chosen, nor any initial population. Only the length- 
scales Sniin/smax and Smax/i are the parameters to be 
decided, where L defines the total volume = L'' of 
an initially empty box. Periodic boundary conditions 
are used so that possible interactions with walls can be 
ignored. In deterministic algorithms, the d + 1 "appro- 
priate" first neighbours must be identified so that the 
center's coordinates of a new sphere could be calculated 
(Soddy's rule); instead this random approach only has 
the non-overlapping constraint, meaning that the overall 
computational complexity is decreased. 
A regular squared mesh with lattice constant aiatticc is 
defined into V and its nodes are used as starting cen- 
ters for the spheres to be inserted. The only requirement 
on the lattice constant is to be sufficiently smaller than 
Smin- A fixed lattice constant could be preferred in some 
cases (the computational cost depends of course on the 
implementation); we tested oiattice = Smin/3 to be a good 
parameter, but a recursive remeshing has been preferred 
and used to assure the requirement to be fulfilled. 
As the first sphere is inserted at random on a starting 
center, it doesn't encounter any other spheres and its di- 
ameter can be expanded to Smax- The starting centers it 
will cover will then be erased. One sphere at a time is 
subsequently inserted, according to the following scheme: 

1: it is placed at random on one of the remaining start- 
ing centers (nodes); 

2: its diameter is increased up to Smax or until it 
touches a previously inserted sphere; 

3: a random displacement within a maximum length 
Ar ^ Smin (average displacement Ar/2) is ac- 
cepted only if this lets the diameter to grow; 

4: step 3 is iterated if the newly calculated diameter 

5: if Smax is reached or if the maximum number of 
displacement attempts (smax/Ar)^ is reached, then 
the procedure stops. 

6: once a sphere has been inserted, the starting centers 
it covers are erased. 

The stopping conditions 5 rely on the possible random 
walk a sphere would need to explore a space of size Smax- 
A maximum displacement length Ar = | Ar| <^ Smin gives 
an average attempted displacement Ar/2 and the num- 
ber Ua of attempted displacements has been chosen to be 
the constant 



(^max \ ^ 
atJ 



4 



corresponding to the number of steps for a random 
walk to explore a region of diameter Smax- Of course a 
higher Ua could result in a higher final packing fraction, 
but our calculation show poor improvement coupled 
with higher computational cost. During the procedure, 
spheres with final s < Smin are erased. The overall 
filling procedure ends when no other starting centers 
are present. In Fig. [2] snapshots are shown for different 



i/Sn 



2 and d — 3. For our simulations 



we used the value L = 1, the others lengths being 
consequently defined with respect to it. Once a value 
is chosen for s,„ax/i, the ratio Smin/smax is investigated. 
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FIG. 2: Snapshots of packings from independent algorithm 
runs: for d = 2 (left) with Smin/Smax ~ 1/5, 1/20, 1/50 for 
A, B and C respectively (and the same Smax/i = 1/4); for 
d = 3 (right) with s^in/smax = 4/10, 2/10, 1/10 for D, E 
and F respectively (and the same Smax/i ~ 1/3). 



2.3 Results 

For the different values of Smax we have tested, we 
observed the same behaviour for decreasing Smin; in the 
next paragraph we report on the results for 370 runs 
performed with Smax = 1/5 in d = 2 and for 170 runs 
with Smax = l/3inc? = 3. These results are shown in 
Fig. [3j where each symbol represents the value averaged 
over 10 independent runs. 



For any configuration of N^^^^^ non-overlapping spheres 
in the total volume V, there always exist some values 
df < d for which the condition ([7| is satisfied. This 
does not imply the existence of a method capable of 
fining volume accordingly to ( [9a] ) and (9b), but simply 
implies that if such an "iterative method" exists, then 
it allows the space to be occupied with a certain df. 
If this method consists in recursively filling the voids, 
each time maximizing the occupied space, then it should 
always present the same fractality as an AP. We expect 
that with such a kind of filling, including the algorithm 
presented here, not only the correct asymptotic power- 
law behaviour has to be obtained, but also that the 



more stringent expressions (9a) and (9b) are fulfilled in 
the overall range of sizes. While power-laws (3a I and 
(3b) are known to work for Smin/smax smaller than 1/5, 



eqs. (9a 



and ( 9b ) are in fact expected to work in the 

max € [0, 



whole interval Smin/smax G [0, 1 
We begin by testing the asymptotic behaviour of our 
osculatory RAP on the power law presented in eq. ^ 
with the values of df for d — 2 and d = 3 obtained 
from previous calculations on AP [T^ [TB]. As can be 
argued from Fig. [3j in both cases good agreement exists 
for low enough values of Smin/smax, proving the correct 
AP asymptotic behaviour for the data obtained by the 
algorithm. 

The numerical results deviate from the fractal asymp- 
totic regime above certain size ratios, where we expect 
instead that our equations should still hold. To this aim 
we fit the data points for d = 2 and d = 3 with cqs. (9a 



and (9b) in the full range of sizes. It is important to 
stress that the value Ng^^^^, which enters the definition 
of N/Ns^^^ and £s,„^^ = 1 — <Ps^^„, is not a fitting 
parameter. In fact it is known as the average number of 
spheres of diameter Smax in the obtained packings. 
Results are shown in Fig. [3j a comparison with the 
values known in the literature is reported in Tabjlj The 
deviation from the asymptotic power-law is evident for 
N/Ns^^^, while the porosity holds its power law form 



(3b) in eq. (9b) 



The iV/A^s,„ax curve results to describe the data better 
then its asymptotic counterpart. The fits finally allow an 
independent estimation of the fractal dimension df: our 
fitted fractal dimensions are in a very good agreement 
with the known values. Note that the possibility of 
fitting in the whole s-range, plus the use of independent 
simulations, allows to evaluate df by simulating a 
relatively small number of spheres. 



Dimension 


df [nmsj 


df (fit) 


d = 2 


1.305684 


1.3045 ± 0.0006 


d = 3 


2.473946 


2.4739 ± 0.0014 



TABLE I: values for the evaluated fractal dimension (in 2 and 
3 Euclidean dimensions) compared with the known ones. 
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based on any a priori size distribution. This algorithm 
aUows to fix, as input parameters, the largest and small- 
est sizes. In the limit of vanishingly small diameters, 
disordered Apollonian packings are recovered with the 
correct asymptotic behaviour, as proved by testing the 
data obtained from the new packing algorithm. The laws 
corrected for finite sizes have been tested by varying the 
interval of sizes for fractal objects (circles in d = 2 and 
spheres and d = 3) and the result of the fits allows to 
verify the values of the fractal dimension which came out 
in agreement with the values known from the literature. 
It is interesting to note that the laws we derived apply to 
the whole range of size ratios. This suggests that even in 
the case of packings with a very narrow interval of sizes, 
the space filling construction preserves its fractal nature. 
These simple results could be applied to the broad class 
of "packing-limited growth" models and physical fractals 
for which the general conditions Q and ([t]) are valid. 
We finally propose to use the rapid convergence to an os- 
culatory packing enhanced with the proposed algorithm, 
together with the possibility of using the whole range of 
sizes for the evaluation of the fractal dimension, as a pos- 
sible feasible test on recent studies on Apollonian gaskets 
at higher dimensionality [26] . 
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FIG. 3: Number of inserted spheres A'' (over Ns„-^^^) as 
function of the scale parameter Smin/smax for d = 2 (top, 
Si-nax/L — 1/5) and d — 3 (bottom, Smin/L — 1/3). Each 
symbol represents the average over 10 independent realisa- 



tions of our random AP. The dashed lines represent eq. ( 3a 
while the solid lines correspond to the fits made using (9a 
Porosities £ for the same syste ms are sh own in the respective 
insets: the curves from eqs. (3b I and (9b I coincide with a 



power law. The positions for samples A, B, C, D, E and F in 
Fig. [2] are indicated on the curves. 



3. Conclusions 

We have studied the properties of space filling packings 
of spheres. In particular we have derived the finite size 
correction to the distribution laws that characterize the 
Apollonian packing fractals. In the case of sizes limited 
to a finite interval, the ratio between the smallest and 
the largest spheres does not go to zero and deviations 
are indeed expected from the typical power laws. Two 
main observables have been studied, the ratio of the to- 
tal number of spheres over the number of largest spheres, 
i.e. N/Ns^^^, and the fraction of unoccupied space, i.e. 
the porosity £. We provide simple analytical expressions 
for them, solely based on the hypothesis that the packed 
spheres totally occupy the space if the minimum diame- 
ter tends to zero. 

In order to test our prediction, an efficient algorithm has 
been introduced to generate osculatory random AP, not 
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